#ifndef _h_models
#define _h_models

#include "geo.h"
#include "tree.h"
#include "model-convex.h"
#include <cassert>
#include <cstdio>

double rand_vx() {
	double vx=(double)(rand()%10000)/10000.;
	//vx=vx*vx;
	vx*=(double)(rand()%2*2-1);
	return vx;
}

struct Light {
	Point p;
	double r;
	Color c;
	Light()
		: r(0) { }
	Light(Point p, double r, Color c)
		: p(p), r(r), c(c) { }
	Point get() const {
		Point cur=p;
		cur.x+=r*rand_vx();
		cur.y+=r*rand_vx();
		return cur;
	}
	Point get(double cx, double cy) const {
		return Point(p.x+r*cx, p.y+r*cy, p.z);
	}
};

struct Model {
	Tree tree;
	vector<Light> lights;
	Model(const vector<Patch> &patches, const vector<Light> &lights)
		: tree(patches) , lights(lights) { }
};

namespace Models {
	
	//void rotate_point_2D(double &x, double &y, double sinA, double cosA) {
	//	double x0=x, y0=y;
	//	x=x0*cosA-y0*sinA;
	//	y=x0*sinA+y0*cosA;
	//}
	
	void rotate_point(Point &p, Point axis, double ang) {
		double x=axis.x, y=axis.y, z=axis.z;
		double s=sin(ang), c=sqrt(1-s*s);
		double R[3][3];
		R[0][0]=x*x+c*(1-x*x), R[0][1]=x*y*(1-c)-z*s, R[0][2]=z*x*(1-c)+y*s;
		R[1][0]=x*y*(1-c)+z*s, R[1][1]=y*y+c*(1-y*y), R[1][2]=y*z*(1-c)-x*s;
		R[2][0]=z*x*(1-c)-y*s, R[2][1]=y*z*(1-c)+x*s, R[2][2]=z*z+c*(1-z*z);
		x=p.x*R[0][0]+p.y*R[0][1]+p.z*R[0][2];
		y=p.x*R[1][0]+p.y*R[1][1]+p.z*R[1][2];
		z=p.x*R[2][0]+p.y*R[2][1]+p.z*R[2][2];
		p=Point(x, y, z);
	}
	
	void rotate(vector<Patch> &patches, vector<Light> &lights, Point axis, double ang) {
		for(size_t i=0; i<patches.size(); ++i)
			for(int j=0; j<3; ++j)
				rotate_point(patches[i].p[j], axis, ang);
		for(size_t i=0; i<lights.size(); ++i)
			rotate_point(lights[i].p, axis, ang);
	}
	
	void rotate(vector<Point> &p, Point axis, double ang) {
		for(size_t i=0; i<p.size(); ++i)
			rotate_point(p[i], axis, ang);
	}
	
	void add_bottom(vector<Patch> &patches, Color cl, double w, double z0) {
		double min_z=INFTY;
		for(size_t i=0; i<patches.size(); ++i)
			for(int j=0; j<3; ++j)
				min_z=min(min_z, patches[i].p[j].z);
		for(size_t i=0; i<patches.size(); ++i)
			for(int j=0; j<3; ++j)
				patches[i].p[j].z+=z0-min_z;
		Point p0=Point(-w, -w, z0);
		Point p1=Point(+w, -w, z0);
		Point p2=Point(+w, +w, z0);
		Point p3=Point(-w, +w, z0);
		patches.push_back(Patch(p0, p1, p2, cl));
		patches.push_back(Patch(p0, p2, p3, cl));
	}
	
	struct Triple {
		int i, j, k;
		Triple(int i, int j, int k)
			: i(i), j(j), k(k) { }
		Triple()
			: i(0), j(0), k(0) { }
	};

	vector<Patch> load_obj(string fname, Color c) {
		FILE *f=fopen(fname.c_str(), "r");
		if(f==NULL) {
			printf("can't open obj file '%s'!\n", fname.c_str());
			return vector<Patch>();
		}
		int n, m;
		if(fscanf(f, "%*s%d%d", &n, &m)!=2) {
			printf("error while reading n, m!\n");
			exit(1);
		}
		// read points and then scale them
		vector<Point> p(n);
		Box box;
		for(int i=0; i<n; ++i) {
			if(fscanf(f, "%*s%lf%lf%lf", &p[i].x, &p[i].y, &p[i].z)!=3) {
				printf("error while reading vertices!\n");
				exit(1);
			}
			swap(p[i].y, p[i].z);
			p[i].y=-p[i].y;
			box.add(p[i]);
		}
		double scale=max(max(box.t.x-box.s.x, box.t.y-box.s.y), box.t.z-box.s.z)/2.;
		if(scale==0.) scale=1.;
		for(int i=0; i<n; ++i) {
			p[i].x=(p[i].x-(box.s.x+box.t.x)/2.)/scale;
			p[i].y=(p[i].y-(box.s.y+box.t.y)/2.)/scale;
			p[i].z=(p[i].z-(box.s.z+box.t.z)/2.)/scale;
		}
		// reading patches and calculating normals
		vector<Triple> tp(m);
		vector<Point> normal(n, Point());
		vector<double> coeff(n, 0.);
		for(int i=0; i<m; ++i) {
			int x, y, z;
			if(fscanf(f, "%*s%d%d%d", &x, &y, &z)!=3) {
				printf("error while reading patches!\n");
				exit(1);
			}
			--x, --y, --z;
			tp[i]=Triple(x, y, z);
			Point cnormal=mult(p[x], p[y], p[z]);
			normal[x]+=cnormal;
			normal[y]+=cnormal;
			normal[z]+=cnormal;
		}
		fclose(f);
		for(int i=0; i<n; ++i)
			normal[i].normalize();
		
		// make patches and return
		vector<Patch> patches(m);
		for(int i=0; i<m; ++i) {
			int x=tp[i].i;
			int y=tp[i].j;
			int z=tp[i].k;
			patches[i]=Patch(p[x], p[y], p[z], normal[x], normal[y], normal[z], c);
		}
		return patches;
	}

	vector<Patch> sphere(Color c) {
		vector<Patch> t;
		const int s_n=5;
		const int s_k=3;
		Point s_p[s_n][s_k<<(s_n-1)];
		double s_h[s_n];
		for(int i=0; i<s_n; ++i)
			s_h[i]=cos(PI/2.*(double)i/(double)(s_n-1));
		for(int i=0; i<s_n; ++i) {
			double y=s_h[i];
			double x=sqrt(1.-y*y);
			for(int j=0; j<(s_k<<i); ++j) {
				double ang=PI*2.*j/(double)(s_k<<i);
				s_p[i][j].x=x*cos(ang);
				s_p[i][j].y=x*sin(ang);
				s_p[i][j].z=y;
			}
		}
		for(int i=1; i<s_n; ++i) {
			for(int j=0; j<(s_k<<i); j+=2) {
				int k0=j>>1;
				int k1=(k0+1)%(s_k<<(i-1));
				int j0=j;
				int j1=(j+1)%(s_k<<i);
				int j2=(j+2)%(s_k<<i);
				#define add(p0,p1,p2) t.push_back(Patch(p0, p1, p2, c)),\
					t.push_back(Patch(p0.neg_z(), p2.neg_z(), p1.neg_z(),c))
				add(s_p[i][j0], s_p[i][j1], s_p[i-1][k0]);
				add(s_p[i][j1], s_p[i][j2], s_p[i-1][k1]);
				add(s_p[i][j1], s_p[i-1][k1], s_p[i-1][k0]);
				#undef add
			}

		}
		return t;
	}
	
	vector<Patch> random_sphere(int n, Color cl) {
		vector<Patch> t;
		Convex3D convex(n);
		srand(123456789);
		for(int i=0; i<n; ++i) {
			double theta=(rand()%10000)/10000.*PI*2.;
			double phi=(rand()%10000)/10000.*PI;
			Point p;
			p.x=cos(phi)*sin(theta);
			p.y=sin(phi)*sin(theta);
			p.z=cos(theta);
			convex.assign(i, p);
		}
		convex.run();
		for(int i=0; i<convex.size(); ++i) {
			Patch patch=convex.get(i);
			patch.c=cl;
			t.push_back(patch);
		}
		return t;
	}
	
	vector<Patch> earth_sphere(int n, int m, Color cl) {
		vector<Patch> t;
		Convex3D convex((2*n-1) * m);
		for(int i=-n+1, cur=0; i<n; ++i)
			for(int j=0; j<m; ++j) {
				double theta=(double)i/(double)n*PI/2.;
				double phi=(double)j/(double)m*2*PI;
				Point p;
				p.x=cos(phi)*cos(theta);
				p.y=sin(phi)*cos(theta);
				p.z=sin(theta);
				convex.assign(cur++, p);
			}
		convex.run();
		for(int i=0; i<convex.size(); ++i) {
			Patch patch=convex.get(i);
			patch.c=cl;
			t.push_back(patch);
		}
		return t;
	}
}

#endif
